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ONE -DIMENSIONAL NUMERICAL ANALYSIS OF THE TRANSIENT 
THERMAL RESPONSE OF MULTILAYER INSULATIVE SYSTEMS 

Claud M. Pittman and Kay L. Brinkley 
Langley Research Center 

SUMMARY 

A one -dimensional numerical analysis of the transient thermal response of multi- 
layer insulative systems has been developed. The analysis can determine the tempera- 
ture distribution through a system consisting of from one to four layers, one of which can 
be an air gap. Concentrated heat sinks at any interface can be included. The computer 
program based on the analysis will determine the thickness of a specified layer that will 
satisfy a temperature limit criterion at any point in the insulative system. The program 
will also automatically calculate the thickness at several points on a vehicle and deter- 
mine total system mass. 


INTRODUCTION 

Some form of analytical tool is necessary for designing thermal protection sys- 
tems (TPS) for high-speed aerospace vehicles. Numerous analyses have been developed 
(refs. 1, 2, 3, and 4, for example); however, these analyses are relatively complex and 
difficult to use and have large computer time and storage requirements. These analy- 
ses were developed to analyze materials systems such as ablative or transpirative heat 
shields. With the advent of projects to develop reusable hypersonic vehicles, such as the 
hypersonic airplane and the space shuttle, interest has shifted to a much simpler reusable 
TPS which does not require such complicated analyses. 

System synthesis computer programs, which determine optimum aerodynamic and 
structural configurations for high-speed vehicles, are also being developed. Because the 
TPS is usually a significant part of the total vehicle mass, an efficient TPS computer pro- 
gram is needed which calculates changes in total heat shield mass with changes in vehi- 
cle configuration and which can be used as a subprogram within the vehicle optimization 
program. 

This paper describes the development and use of a one- dimensional, multilayer 
insulative TPS computer program. Both the differential and finite-difference equations 
are given. The computer program input quantities are listed and the computer program 



output for a sample case is presented. The surface environment, initial temperature 
distribution, and material properties must be specified. 

SYMBOLS 


*w 




X 




Any consistent set of units can be used in this analysis. 

heat capacity of heat sink, fpcnt) , j/m^-K 

' ^ 'heat sink 

specific heat, J/kg-K 


total enthalpy at edge of boundary layer, J/kg 

local enthalpy of fluid at front surface temperature, J/kg 

effective heat transfer coefficient, W/m2-K 


effective heat transfer coefficient to fluid at back surface, W/m2-K 

effective heat transfer coefficient at air gap, W/m2-K 

thermal conductivity, W/m-K 

cold- wall convective heating rate, W/m2 

radiant heating rate, W/m2 

temperature, K 

temperature to which back surface radiates, K 
temperature of fluid at back surface, K 
heat sink thickness, m 
coordinate normal to surface 
thickness of jth layer, m 
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a absorptance of front surface 

e effective emittance 

effective emittance of back surface 
€g effective emittance at air gap 

6g effective emittance of front surface 

p density, kg/m^ 

0 Stefan-Boltzmann constant, W/m^-K^ 

T time, sec 

Subscripts: 

1,2,3 first, second, and third stations or layers 

b back surface 

g air gap 

1 interface 

j layer 

L last station 

n integer 

s front surface 

In the finite-difference equations, primed quantities are evaluated at the end of the 
time step. Double primed quantities are functions of temperature but are evaluated at the 
end of the time step and therefore require iteration. 
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GENERAL DESCRIPTION 


This analysis was developed primarily to analyze radiative heat shield systems. It 
can be used for most systems in which the material surfaces do not move with time and 
in which two-dimensional effects can be neglected. The sailient features of the analysis 
are as follows; 

(1) From one to four layers of different materials can be considered. More lay- 
ers can be included by simply changing the dimension statements in the first common 
block, but this change will increase computer storage requirements. 

(2) One layer can be an air gap, but the air gap must have a layer on both sides of 
it. Heat transfer across the air gap can be by radiation, convection, and conduction. 

(3) Perfect thermal contact between layers is assumed. Contact resistance can be 
simulated by using a separate layer or the air gap option at interfaces which have signifi- 
cant contact resistance. 

(4) In each layer, the density is constant, the specific heat can be a function of 
temperature, and the thermal conductivity can be a function of both temperature and 
pressure. 

(5) If any layer has a thermal conductivity much larger than the conductivity of the 
other layers, this layer can be input as a concentrated heat sink. Concentrated heat sinks 
absorb heat according to their material properties and thickness (pcp^ but do not have a 
thermal gradient nor occupy space. Concentrated heat sinks can be used at the front sur- 
face, at each interface, and at the back surface. 

(6) Both convective and radiant heat inputs to the front surface can be used. If 
cold-wall convective heating is used, a hot-wall correction (based on enthalpy) can be 
used. If radiant heating is used, the surface absorptance can be specified as a function 
of temperature. 

(7) Radiation from the front surface, across an air gap, and from the back surface 
can be used. Radiation from the back surface is to a specified temperature which can be 
a function of time. All emittances can be functions of temperature. 

(8) Heat transfer by convection can be used at the back surface. The heat transfer 
coefficient can be a function of the back surface temperature and the temperature of the 
fluid can be a function of time. 

(9) The temperature at any one location can be specified as a function of time. 

(10) The computer program has been extended to permit the calculation of the total 
vehicle heat shield mass with one set of input data. Details of this option are given in a 
following section. In brief, the program calculates heat shield mass corresponding to 
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specified temperature limits at as many as 10 points around the vehicle and integrates 
over these points to obtain an average heat shield mass per unit area for the entire 
vehicle. 


ANALYSIS 

The equations describing the transient temperature response of a one-dimensional, 
multilayer, insulation system have been developed, put in finite -difference form, and pro- 
gramed for solution on the digital computer. The coordinate system originates at the 
front surface. 


Governing Differential Equation 
The governing differential equation for any layer is 




( 1 ) 


The density p within a layer is constant, the specific heat Cp can vary with tempera- 
ture, and the conductivity k can vary with temperature and pressure. The initial tem- 
perature, as a function of position, must be specified. 


Boundary Conditions 

Energy balances are used at each boundary. 

Front surface .- The front surface boundary condition is 


Cold-wall convective Hot-wall Radiative Conduction Rerad at' Heat absorbed by 

heating rate correction heat input to interior ^ ^ surface heat sink 



The two quantities in the hot-wall correction (the enthalpy of the fluid at the front surface 
temperature hw and the enthalpy at the edge of the boundary layer he) should be based 
on absolute zero. If the base is room temperature, numerical instabilities can result if 
hg approaches zero. 

Interface . - The boundary condition at the interface between two layers of material 
is 


J ax 




( 3 ) 


5 



Back surface.- The back surface boundary condition is 



where is the effective emittance of the back surface and hf is the effective heat 
transfer coefficient at the back surface. 

Air gap.- The boundary condition at the front side of the air gap is 


-k.^| 

J 8x 




_-k II 

S 8x 




+ o^eglT 




+ h 


Ti - T^ 


i+1 




8Ti 

'dT~ 


(5) 


* 

where e„ and hg. are effective values of the emittance at the air gap and the heat trans- 

& o 

fer coefficient at the air gap, respectively. 

At the back of the air gap, the boundary condition is 




+ aeg 








+ 



( 6 ) 


Finite -Difference Equations 

The finite-difference stations are spaced at equal intervals throughout a given layer, 
as shown in figure 1. The first layer includes a station at the front and the back of the 
layer. All other layers include a station at the back of the layer. The air gap has no 
stations within the layer. 

The differential equations are put in finite-difference form by the use of Taylor's 
series expansions. Forward, central, and backward differences are used. At the bound- 
aries, the second derivative in the governing differential equation (eq. (1)) is partially 
expanded to first derivatives and the boundary condition is inserted in place of the first 
derivative at the boundary. Thus, both the governing equation and the boundary condition 
are satisfied at the boundary. In the equations, all terms which are not explicit functions 
of time are averaged over the time step because a Taylor's series expansion of the time 
derivative at At/2 shows that this procedure makes the solution correct to At 2. (See 
ref. 5.) 

A summary of the finite-difference equations is given in appendix A. These equa- 
tions yield an essentially tridiagonal matrix of unknown temperatures. A procedure 
based on Gauss' elimination method is used to solve the matrix. 
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Heat Shield Mass Calculation 

The computer program has an option which permits iteration of the thickness of a 
specified layer to satisfy a specified temperature limit criterion. For example, the min- 
imum thickness of the primary insulation layer, which will limit the back surface to a 
specified temperature, can be determined. Any one layer can be specified for the thick- 
ness iteration and any temperature limit can be imposed on any one station. 

The program also has an option to determine minimum thicknesses at as many as 
10 different locations, which are equivalent to different body points around the vehicle. 

The surface environments (cold-wall convective heating rate and pressure) at body points 
2 to 10 are specified as fractions of the initial cold-wall convective heating rate and pres- 
sure tables. These fractions can be functions of time to accommodate transitions to tur- 
bulent flow. The enthalpy is assumed to be constant around the body. 

The heat shield mass at each body point can also be calculated. In this calculation, 
the mass of each layer and a specified mass (through which the mass of concentrated heat 
sinks can be added) are included at each body point. 

By assuming that the surface environment at each body point corresponds to an 
average environment over a specified percentage of the total vehicle area, an average 
heat shield mass per unit area for an entire vehicle can be computed. 

COMPUTER PROGRAM 

This computer program (Langley designation A4596)^ was written in FORTRAN 
language for the Control Data 6000 series digital computer under the SCOPE 3.0 operating 
system. The program requires approximately 45000 octal locations of core storage. A 
flow chart for the program is shown in figure 2. The program input quantities are listed 
and defined in appendix B. 


Computer Program Accuracy 

The accuracy of the computer program was determined by comparing the results 
from the numerical solution with two exact solutions and also comparing results with the 
established computer program of reference 1. 

The first exact solution used was taken from reference 6. Equation (All) of refer- 
ence 6 gives the transient temperature response of a constant property slab, subjected to 
a constant heating rate on one side and perfectly insulated on the other side. For the 
case selected, the average surface temperature rise rate was 740 K/sec. The maximum 


^This program is available as LAR- 12057 from COSMIC, Computer Center Informa- 
tion Services, 112 Barrow Hall, University of Georgia, Athens, Georgia 30601. 
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errors for several numerical solutions are given in table I. In all cases the maximum 
error occurred at the first printout (0.25 sec). (The time step was 0.03125 sec.) The 
accuracy of the solution was not significantly affected by the number of layers used. The 
smaller error criterion increased accuracy by a factor of 100 but increased computer 
time by a factor of about 10. 

The exact solution was also computed with the air gap included to check the air-gap 
equations for correctness and accuracy. When the air-gap equations were used, the error 
increased slightly and decreasing the error criterion did not increase the accuracy sig- 
nificamtly. The increased error was probably caused by the less accurate two station 
forward and backward differences used at the air gap. 

The other exact solution was taken from reference 7. This exact solution was for a 
constant property slab with a constant surface temperature, losing heat at the back sur- 
face to a constant-temperature fluid. For this case, the maximum error in the steady 
state solution was 0.002 percent. 

The computer program was compared with the analysis of reference 1 by using a 
variable property, transient case. Both programs gave essentially the same results with 
the temperatures differing in about the eighth place. 

Sample Case 

The computer program output for a sample case is given in appendix C. The case 
was selected to demonstrate, fairly simply, the versatility of the program. The environ- 
ment roughly represented a flight test of a hypersonic research airplane. The heat shield 
consisted of a metallic outer skin over a metal -foil -encapsulated fibrous insulation with 
an air gap between the inner foil surface and the main structural skin. This heat shield 
system was modeled for the computer program as shown in figure 3. A temperature limit 
criterion for the main structural skin was specified and the insulation layer thickness was 
iterated until this criterion was satisfied. Calculations were made at 6 body points around 
the vehicle at nominal heating rate ratios (computer input quantity QRATT) of 1.0, 0.8, 0-6, 
0.4, 0.2, and 0.1. The nominal QRATT values of 0.4, 0.2, and 0.1 were increased during 
the latter part of the heat pulse to simulate transition to turbulent flow. The heat shield 
unit mass was integrated around the vehicle as a function of percent vehicle area at each 
heating rate ratio. The total heat shield mass can be found by multiplying the final unit 
mass value by the total vehicle area. 

A detailed description of the sample case printout is given in appendix C. The 
results from the sample case are given in figures 4, 5, 6, and 7. Figure 4 shows the 
insulation layer thickness iterations for the specified nominal QRATT values. For the 
initial QRATT value of 1.0, 3 iterations were required to satisfy the back surface tem- 
perature limit criterion, mainly because of the inaccurate initial thickness estimate. The 
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extrapolation equation used to estimate thickness at the lower QRATT values generally 
gave good estimates. Because the extrapolation equation is based on the nominal QRATT 
values, the equation tended to underestimate the thickness requirements of the QRATT 
values (0.4, 0.2, and 0.1) which were increased during the later part of the heat pulse to 
simulate transition to turbulent flow. The last calculated insulation layer thickness at 
each QRATT value is plotted in figure 5 as a function of the nominal QRATT values. 

The assumed heating rate distribution over the vehicle is given in figure 6, which 
shows, for example, that 10 percent of the vehicle area has an average heating rate equal 
to the initial convective heating rate (QRATT = 1.0) and 15 percent of the vehicle has an 
average heating rate equal to 80 percent of the initial heating rate. Figure 7 shows the 
cumulative average heat shield unit mass around the vehicle as a function of percent vehi- 
cle area. Each mass point includes the unit mass of the layers and the unit mass of the 
concentrated heat sinks. The average unit mass for the vehicle is 7.7 kg/m^. This unit 
mass multiplied by the total vehicle area gives the total heat shield mass. 

CONCLUDING REMARKS 

A numerical analysis of the transient response of multilayer insulative heat shield 
systems has been developed. The differential equations, finite-difference equations, 
computer program input listing, and computer program output for a sample case are 
given. 

The computer program is versatile and efficient and has relatively small computer 
storage requirements. Within a single calculation (for a specified trajectory and vehicle 
heating distribution), the program will iteratively determine the heat shield thickness 
required for a specified temperature limit criterion, determine the thickness require- 
ments at a number of points on the vehicle, and determine the average heat shield unit 
mass for the vehicle. The program is suitable for use as a subroutine in large computer 
programs for vehicle system synthesis. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
March 11, 1976 
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APPENDIX A 


FINITE-DIFFERENCE EQUATIONS 

The differential equations were put in finite-difference form by the use of Taylor's 
series expansions. Three station forward, central, and backward differences were used. 
The terms derived from the second derivative of equation (1) were averaged over the time 
step making the solution correct to order At2. 

At the boundaries, the second derivative in the governing differential equation 
(eq. (1)) was partially expanded to first derivatives and the boundary condition was in- 
serted in place of the first derivative at the boundary. Thus, both the governing equation 
and the boundary condition are satisfied at the boundary. 

The finite-difference form of equation (1) and combinations of equation (1) and the 
boundary conditions (eqs. (2) to (6)) are presented in this appendix. In the following equa- 
tions, unprimed quantities are either evaluated at the beginning of the time step or are 
constant with time. Primed quantities are evaluated at the end of the time step. Double 
primed quantities are functions of temperature but are evaluated at the end of the time 
step and therefore require iteration. 

The first subscript on material property values refers to the layer. Further sub- 
scripting indicates the station temperatures at which the material property is evaluated. 
Thermal conductivities are evaluated at the average temperature between two stations. 


Layer Equation 

The finite -difference equation for any layer is 



2 Axj2 2 Axj2 



(Al) 
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APPENDIX A 


First Station 

The finite-difference equation at the heated surface is 



Interface Equation 

The finite -difference equation at any solid-solid interface is 




■i-l 


2 AXj 




T, 


2 Ax. 


i+i 


- "i-2) - - ^1-2 

6 Ax^ 


4 C,(T. - Ti 


i+l,i+2 


6 Axj 


i+1 


3 At 


p. AXj 


(‘=P.i)l 


4 + C„ , 




2 At 




'P,3+l)i ^ (^P,3+l)i 


T! - T. 
1 1 


2 At 


= 0 (A3) 
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Back Surface 

The finite-difference equation at the back surface is 


5 ^L,L-1 




2 Ax4^ 


L 1 .L- 2 V ^ - ae;;(Tn4 

6 Axj 2 3 Ax^[ b( Lj L bl bj 


+ ae,,fT. 4 - tA + (ht) hi - A + hthj - T,'\ I - 


^b\^L - ^b^-^^f; ^fj-“f^^L ^f) 3AX4 at 


("pjIl rpdjLjrL - Tl 


Air Gap 

The finite -difference equation at the front side of the air gap is 


SAXj 


t: 1 - T o + k, 


■3/i-ii-2r-i ^-2rn;i_ii-2ri-i ^- 2 , 


T, . - T 4 


8 Axj 


T;.i-T;u(kg) fT.^i-T.) 


- ^ T: - T[h . hg^i T. - T.^i) 


“g,i’(Tl’f t; - + <JejT* - tQ . 2C, 


3 Ax.p. I'P.ljl " 5.i), Ti ■ 


^ ^ ’'j+lPj+1 
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At the back side of the air gap, the finite-difference equation is 






Ax 


J 


8 Ax, 


3+1 


k 


Tl.o - tI.A + k, 




8 Axj+i 


■-y m-i - t; - - T, - 


r„T! \3„t 

Ti-1 T,_i 


T”fT: 


"'g.ir 


T - Tj 

+ 2 Ci 


^pjJi ^ rpd/1 


t:-t. 


2 At 




‘^p,j+i)i ^ I'^p^i+iji 


T. -Ti 


= 0 
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APPENDIX B 


COMPUTER PROGRAM INPUTS 

The input data quantities for the computer program are defined in this appendix. 

All data storage locations are zeroed out before the data are used. Thus, if any quantity 
is zero for a specific case, the quantity can simply be left out of the data input for that 
case. 

Input quantities which are functions of time or temperature are used in tables. 

Each table has four quantities. As an example, consider the table for enthalpy hg. 

(1) The quantity MHE is the order of interpolation. When MHE = 0, hg is 
constant. When MHE = 1, linear interpolation is used. When MHE = 2, second- 
order interpolation is used. 

(2) The quantity NHE is the number of values in the hg table. When NHE = 0, 
hg is constant. 

(3) TAUHE lists the time values corresponding to each value in the hg table. 

When TAUHE = 0, hg is constant. This independent variable table must be either 
monotonically increasing or decreasing. 

(4) HET lists the hg values. 

From this discussion it follows that when hg is constant, only the quantity HET need be 
specified. 

HET must always have a nonzero value because it appears in the denominator in 
equation (2). If the QC table is being used, QRATT must have a value because it multi- 
plies QC. 

After each input quantity which can have more than one value, the maximum number 
of values is given parenthetically. When the letter N appears in the parentheses, either 
alone or as a second dimension, N refers to the layer in question and cannot exceed 4. 

In the following listing, whenever possible, the symbol used in equations (1) to (6) is 
given after the computer symbol definition. The listing is arranged semialphabetically. 
The quantities describing a table (usually four) are grouped together and the name of the 
dependent variable fixes the alphabetical location of the table. 

Front Surface Absorptance 
MAL order of interpolation for ALT 

NAL number of values in ALT 
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TTAL(IO) 

ALT(IO) 

AREA(IO) 

CF 

CIT(N,10) 


MCP(N) 

NCP(N) 

TTCP(10,N) 

CPT(10,N) 

CST(IO) 

DTAUl 

ENDTAU 

ME PS 
NEPS 
TTEPS(IO) 
EPST(IO) 

MEPSB 

NEPSB 

TTEPSB(IO) 

EPSBT(IO) 


temperature table for ALT 
absorptance of front surface a 


percent of total vehicle area to be included in each body point heat 
shield mass calculation 


convergence factor for iterating on layer thickness (if convergence is 
too slow, make CF > 1.0; if thickness iteration diverges, make 
CF < 1.0; if CF is not specified, CF will equal 1.0) 


(pcpt) 

point 


of heat sinks at interfaces Cj 


(one set of values for each body 


Layer Specific Heats 
order of interpolation 
number of values 
temperature table 
specific heat Cp 

^pcptj of heat sink at front surface Cg (one value for each body point) 

initial time step At (the time step is increased or decreased auto- 
matically during the calculation depending on the number of iterations 
required for convergence) 

time at which computation stops 

Front Surface Emittance 
order of interpolation 
number of values 
temperature table 
emittance of front surface €g 

Back Surface Emittance 
order of interpolation 
number of values 
temperature table 
emittance of back surface ej^ 
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MEPSG 

NEPSG 

TTEPSG(IO) 

EPSGT(IO) 

ERR0R 


MHE 

NHE 

TAUHE(20) 

HET(20) 


MHF 

NHF 

TTHF(IO) 

HFT(IO) 


MH 

NH 

TTH(IO) 

HT(IO) 


MHW 

NHW 

TTHW(15) 

HWT(15) 


Air-Gap Emittance 
order of interpolation 
number of values 
temperature table 
emittance at air gap 6g 

acceptable error in temperature convergence, (T' - T")/T* (0.001 is 

usually acceptable) 

Total Enthalpy at Edge of Boundary Layer 
order of interpolation 
number of values 
time table 

total enthalpy at edge of boundary layer hg 

Back Surface Heat Transfer Coefficient 
order of interpolation 
number of values 
temperature table 

heat transfer coefficient to back surface hj 

Air-Gap Heat Transfer Coefficient 
order of interpolation 
number of values 
temperature table 

heat transfer coefficient at air gap h* 

Enthalpy of Fluid at Front Surface Temperature 
order of interpolation 
number of values 
temperature table 

local enthalpy of fluid at front surface temperature h^ 
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I(N) 

number of stations in each layer (first layer must have at least 
3 stations; other layers must have at least 2 stations; air gap has 
only 1 station) 

Air gap 

1ST AS 

number of station at back surface of air gap, 1 + ) I(N) 

n=l 


Thermal Conductivity 

NP(N) 

number of pressure values in each k table (if NP(N) = 1, 
k = f(T) only) 

TABP(5,N) 

pressure table 

TABT(10,N) 

temperature table 

NK(N) 

number of k values in each table (NK/NP = Number of temperature 
values in each table) 

TABK(50,N) 

thermal conductivity values k 

LCHAN 

layer whose thickness is changed to satisfy a temperature limit 
criterion 

NCHAN 

station at which temperature limit criterion is specified 

TCHAN 

temperature limit at station NCHAN which will be satisfied by 
iterating on thickness of layer LCHAN 

NLYERS 

number of layers 

NQP 

number of body points (points around vehicle) to be analyzed (equal 
to number of PRATT and number of QRATT values) 

PFREQ 

print frequency 


Body Point Pressure Ratios 

MPRAT(IO) 

order of interpolation 

NPRAT(IO) 

number of values 

TAUPRAT(10,10) 

time-body point table 

PRATT(10,10) 

fraction of PT values, used to obtain pressure dependence of 


thermal conductivity k at different body points 
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Pressure 

MPT order of interpolation 

NPT number of values 

TAUP(IO) time table 

PT(IO) pressure history used to calculate conductivity values 


Cold- Wall Convective Heating Rate 
MQC order of interpolation 

NQC number of values 

TAUQC(20) time table 

QCT(20) cold-wall convective heating rate 


Body Point Heating Rate Ratios 
MQRAT(IO) order of interpolation 

NQRAT(IO) number of values 

TAUQRAT(10,10) time -body point table 

QRATT(10,10) fraction of QCT values as a function of time and body point 


MQR 

NQR 

TAUQR(20) 

QRT(20) 

RH0(N) 

SIGMA 

T(IOO) 

TAU0 


Radiant Heating Rate 
order of interpolation 
number of values 
time table 

radiant heating rate qj^ 
density p 

Stefan- Boltzmann constant a 
initial temperature distribution T 
starting time 


Temperature to Which Back Surface Radiates 
MTB order of interpolation 

NTB number of values 
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TAUTB(IO) 

TBT(IO) 

MTF 

NTF 

TAUTF(IO) 

TFT(IO) 

NTSP 

MTS 

NTS 

TAUTS(20) 

TST(20) 

X(N) 

WTCl(lO) 


time table 

temperature to which back surface radiates T^^ 

Back Surface Fluid Temperature 
order of interpolation 
number of values 
time table 

temperature of fluid at back surface Tf 
station at which temperature is specified 


Station Temperature History Specification 
order of interpolation 
number of values 
time table 

specified temperature history at station NTSP 


V 


layer thickness 

total heat sink unit mass, ^ pt, to be added to layer unit masses to 
get total heat shield unit mass at each body point (one value for each 
body point) 
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SAMPLE CASE 

The computer printout for a sample case is given in this appendix. The first infor- 
mation given is a list of the input data. The entire dimensioned field of each input quan- 
tity is given. Because the entire field for each quantity is seldom filled, a great many 
zeros appear in the input listing. The order of the input list is approximately as follows: 
program control data, material property data, environment data, and data for calculations 
at several body points. 

Following the data input list, the regular printouts begin. The first printout occurs 
after the solution has converged for the first time step. Subsequent printouts are at the 
specified print frequency. For the sample case, the print frequency was made much 
larger than usual to reduce the number of pages of printout. The computer always prints 
out the last calculation for any normal program stop. For each thickness iteration in 
this sample case, the last printout is at the time when the back surface temperature 
starts to decrease. 

Each printout starts with the time and the time step. Next, the temperature at each 
station (starting at the front surface) is given. Each line contains six temperatures. 
Therefore, if a layer contains more than six stations, the temperatures cover more than 
one line. However, the temperatures for the next layer always start on a new line. The 
last temperature in each layer and the first temperature in the next layer are identical 
because this station is common to both layers. 

The rest of the quantities in a printout are quantities which vary with time, thick- 
ness iteration, or body point. These quantities are defined as follows: 

QC cold-wall convective heating rate at the beginning of the time step 

QCP cold- wall convective heating rate at the end of the time step 

QR radiant heating rate at the beginning of the time step 

QRP radiant heating rate at the end of the time step 

QRR reradiation from the front surface aCg at the beginning of the 

time step 

HE total enthalpy at the beginning of the time step 

HEP total enthalpy at the end of the time step 

HWO enthalpy of the fluid at the front surface temperature (at the beginning 

of the time step) 
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HW 

X 

QAERO 

QAEROP 

TF 

QRAT 

PRAT 

WTCI 

AREA 

WTBP 

WTSUM 


enthalpy of the fluid at the front surface temperature (at the end of the 
time step) 

thickness of the layer whose thickness is being iterated 

net heat input to the front surface at the beginning of the time step 

qc(l - hw/he) + 

net heat input to the front surface at the end of the time step 

- ^V^e) + 

temperature of the fluid at the back surface (at the beginning of the 
time step) 

heating rate ratio at the end of the time step 
pressure ratio at the end of the time step 
total heat sink unit mass at this body point 

fraction of the total area which has an average heating environment 
corresponding to this body point 

sum of the layer unit masses at this body point (this value does not 
change until the layer thickness iteration has converged) 

cumulative average unit mass / (WTCI + WTBP)AREA (this value 

Lj 

does not change until the layer thickness iteration has converged) 


This sample case was run with 45000 storage locations and the execution time was 
120 sec on a CDC 6400 computer. 




Input Dat 

$DATAl 



TAUn 


o 

• 

o 

OTAUl 

= 

0.3125E-01, 

ENOT AU 


0. le+OA , 

PFRFO 

= 

0.2'= + 03, 

ERROR 


0. lE-02 , 

SIGMA 

= 

0, S6696 LE-IO, 

NLYERS 



I 


20, 1, 2, 0, 

X 


0.30A8E-01, 0*3D^8E-02, D.509E-03 

ISTAS 


21, 
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RHO 

is: 

0.96F+02, 

0.3E + 00, 0. 

27231E+04, 

o 

. 

o 




MCP 

= 

1, 

1, 1, 

0 , 








NCP 

= 

2. 

1, 1, 

0, 








TTCP 

- 

0.2T5F+03, 

0.167E+04, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 0.0 



0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 0.0 



0*0f 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 0.0 



0. Of 

O.Of 

0.0, 

0.0, 

0.0, 

0.0,, 

0.0, 

0.0, 

O.Of 


CPT 

— 

0.62 

■^6E4-00, 

0.12552e+01t 0. 

0, 0. 

0, 0. 

0, 0. 

0, 0. 

0, 0.0, ( 



O.Ot 

0.1046E+OI, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 

0.0, 

0.0, 0.1 



O.Of 

0.9205E+00, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, o.< 



O.Of 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 

0.0, 0.0 

MAL 


Of 










NAL 

= 

0, 










TTAL 

= 

0.0, 

0.0, 

0.0, 

O.Of 

O.Of 

o 

• 

o 

0.0, 

0.0, 

0.0, 

0.0, 

ALT 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

o 

« 

o 

0.0, 

0.0, 

MEPS 


0, 










NEPS 


0, 










TTEPS 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

EPST 

= 

0.8E+00, 0 

.0, 0 

.0, 0 

.0, 0 

.0, 0 

. 0, 0 

.0, 0 

.0, 0 

.0, 0.0, 

MH 


0, 










NH 


0, 










TTH 

= 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

HT 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

MEPSG 


0, 










NEPSG 


0, 










TTEPSG 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 

0.0, 

EPSGT 


0.8f: + 00, 0 

. 0, 0 

.0, 0 

o 

o 

. 

.0, 0 

.0, 0 

.0, 0 

. 

o 

o 

. 

o 

o 

. 

o 

MHF 


0, 










MHF 


0, 










TTHF 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

o 

• 

o 

0.0, 

0.0, 

HFT 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

MEPSB 


0, 










NEPSR 


0, 










TTEPSB 

= 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

o 

. 

o 

0.0, 

0.0, 

EPS8T 

= 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

NP 


3, 

1 , 1, 

0, 
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TABP 


o«ie-oi, 0 

. lE+00, 

0. lE^-01 , 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 



0«0? OaOf 

0.0, 

0.0, 

0.0, 

0.0, 

OoO, 

0.0, 

0.0, 

0.0, 

0.0, 

TABT 


0-275H+03, 

0. 167E + 04, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 

0.0, 

0.0, 



0« 0 t 0 • 0 1 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



OoO» o«o» 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 

0.0, 

0.0, 



0*0» 0«0» 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



NK 

= 

6, 1, 1, 

0, 









TASK 

= 

0,312E-04, 

0.2496E-03, 

0.4368E-04 

, 0.3744E-03, 0 

.49928- 

04, 



0.^368E-03, 

0.0, 

0.0, 

0.0 

f 0.0, 

0.0, 

0.0, 

0. 0 

, 0.0, 

o 

. 

o 



0*0, O.Of 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0«0» 0«0f 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 



0*0» 0®0f 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0«0t 0«0» 

0.312E 

-04, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0» 0®0» 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0«0» 0«0» 

0.0, 

0.0# 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0®0, O.Ot 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 

0.0, 

0.0, 



OoOy 0»0j 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.230888+00 



0®0» 0»0f 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0«0» 0«0f 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0®0t 0«0t 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0®0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



O.Of 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 

0. 0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

MPT 

= 

1, 










NPT 

= 

4, 










TAUP 

ss 

0.0, 0.7E+02, 0. 

22E+03, 0 

.lE+04, 

0.0, 

o 

. 

o 

0.0 

o 

. 

c 

o 

. 

o 



0.0, 










PT 


0.6E-01, 0 

.25E-01 

, 0. 

lE+00 

, O.lEfOl, 

0.0, 

0.0, 

o 

• 

o 

0.0, 



0.0, 0.0, 










T 

= 

0.294E4-03, 

0.294E+03, 

0.294E+03 , 

0.294E4-03, 

0.294E^'03, 




0.294E+03, 

0.294E4-03, 

0.294E + 03 , 

0.294E+03, 

0. 294Ei*03, 




0.29^^E4•03, 

0. 294E+03, 

0.294E+03 , 

0.294E^‘03, 

0.294E+03, 




0.294E4-03, 

0. 294F4-03, 

0.294E+03 , 

0.294E+03, 

0.2948^-03, 




0.294Ef03, 

0.294E+03, 

0.294E+03, 

0.294E-I-03, 

0.294E+03, 




0.294E^03, 

0.294E+03, 

0.294E+0? , 

0.294E4-03, 

0.294EH-03, 




0.2^4E<-03, 

0.294E+03, 

0.294F+03, 

0.294E+03, 

0.294E+03, 




0.294E^'03, 

0.294E+03, 

0.294E+03, 

0.294E^-03, 

0.2948^-03, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0. 0, 

0.0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0, 0.0, 

0 . (3 , 

0.0, 

O.Of 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 



0.0, 0.0, 

0.0, 

0.0, 








MHW 


1, 










NHW 


10, 










TTHW 


0.278E1-03, 

0. 555E + 03 , 

0.833E+03 , 

0. llllE+04, 

0. 13898+04 

, 



0.1667E+04, 

0. 1944E+04 

, 0. 

2222E^04, 0. 

25E+04 

, 0 « 

27788 + 04, 0. 



0. 0 , 0.0, 

0.0, 

0.0, 
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HWT = 0.267F + 03, 0.548F + 07, 0.844E+03. 0. ll6rL<-04, 0. 1489F + 04i 

0.lS33F+04» 0.2188F + 04i 0. 2564'^+0‘tj 0.3011F<-04» 0.3269F + 04» 

0*0f OaOt 0*0t 0*0f 0»0f 

MQC = 1 , 

NQC = 15, 


TAIJQC = 0.0, 0.26E1-02 , 0.45E + 02, 0.6lF.<-02 , 0. 66F. + 02, 0.74F + 02, 

0.88F4-02, 0.92E + 02, 0.103E + 03, 0.114E+03, 0. 142E + 03, 0.17F + 03, 
0.197E4-03, 0. 229E + 03, O.lF+04, 0.0, 0.0, 0.0, 0.0, 0.0, 

QCT = 0.1135E+02, 0.3405E+02, 0.7945F+02, 0.1816F+03, 0.26105F+03, 

0.3405F+03, 0.5675F+03, 0.5675E+03, 0.454Ff03, 0.3405E+03, 

0.227E+03, 0. 1135E + J3, 0. 227E+02, 0.U35F + 02, 0.454Ff01, 0.0, 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 

MQR = 0, 


NQR = 0 , 


TAUQR 

o o 

9 9 

o o 

M 

0.0, 

0.0, 

o o 

• • 

o o 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

o o 

» « 

o o 

0.0, 

0.0, 

0.0, 

o 

• 

o 

QRT 

- 0*0t 

0,0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

o o 

* • 

o o 

0.0, 

0.0 


MHE- = 1 , 


NHE =15, 

TAUHE = 0.0, 0.26E + 02, 0.45F>02, 0.61E+02, 0.66E4-02, 0.74E+02, 

0.88C«-02, 0.92F-t-02, 0.103E + 03, 0.114E + 03, 0. 142E + 03, 0.17E + 03, 
0.l97E-t-03 , 0.229F + 03, O.lE+04, 0.0, 0.0, 0.0, 0.0, 0.0, 

HET = 0.397E+03, 0. 629E<-03, 0.954E + 03, 0. 1366E + 04, 0. 158E-t-04, 

0.1868E+04, 0.246E+04, 0. 246E+04, 0.2153E+04, 0.1868E+04, 

0.1357E+04, 0.937E + 03, 0.617F+03, 0.395E+03, 0.29E4-03, 0.0, 




0.0, 

0.0, 

o 

« 

o 

o 

. 

o 







MTB 

= 

0, 










NTB 


0, 










TAUTB 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

TBT 

= 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

MTF 


0, 










NTF 


0, 










TAUTF 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

TFT 

= 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

NTSP 

= 

0, 










MTS 


0, 










NTS 


0, 










TAUTS 


0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 

0.0, 
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TST 

■s: 

O.Oi O.Ot 

O.Oi 0.0, 

0.0, 0.0, 

0.0, 0.0, 

TCHAN 

= 

0.422E*h03, 


LCHAN 


If 


NCHAN 

= 

23, 


NQP 

= 

6, 


CF 


O.C, 


NQRAT 

= 

1, 1, If 

4 , 5 , 6 , 

MORAT 


o 

c 

o 

1, 1, l» 

TAUQRAT 


0.0, 0.0, 

0.0, 0.0, 

0.0, 0.0, 

0.0, 0.0, 


0 . 0 , 

0 . 0 , 


0 . 0 , 

0 . 0 , 


0 . 0 , 

0 « 0 f 


O.Of 
0 « 0 1 


0«0, 

0.0« 


0-0, 


0«0, 


Ot Oi 


QRATT 


MPRAT 

NPRAT 

TAUPRAT 


PRATT 


0 . 0 , 0 . 0 , 
0. 12F+03, 
0. 103^:^*03, 


0 . 0 , 


0, 

0 . 0 , 

0 . 0 , 

0 . 0 , 


0, 


0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 


0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 


0 . 0 , 

0.1F4-0A, 0.0, 0.0, 

0. llF+03, 0. 12E+03, 

0.P2E-«’02, 


0 . 0 , 0 . 0 , 
0 . 0 , 0 . 0 , 
0 . 0 , 0 . 0 , 
0 . 0 , 0.0 
O.lE+04, 


0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 
0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 
0.0, 0.0, O.llAF+03, 

0 . 0 , 0 . 0 , 0 . 0 , 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 


0 . 0 , 

0 . 0 , 

0 . 0 , 

O.Qi 


0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 


O-IE^-OB, O.llF + O?, 0.1BE^03, O.IF^-OA, 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 


O.Of 

0 . 0 , 

0 . 0 , 


0. lE+Ol , 
0.8E+00, 
0.6F^-00, 
O.^F+00, 
0 . 0 , 0 • 
0 . 0 , 0 . 
0.15F+00 
0 . 0 , 0 . 
0 . 0 , 0 . 
0 * 0 , 0 . 
0 . 0 , 0 • 

0 , 0 , 


0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 

O.AE + 00, 0.5E4*00, 0.4E+00, 

2E + 00, 0.2E4-00, 0.25E4'00, O 

0, 0.0, O.IE^-OO, 0.1E4'00, 


0.0, 

0.0, 

0.0, 

0.0 

0.0, 

0.0, 

0.0, 

0.0 

0.0, 

0.0, 

0.0, 

0.0 

0.0, 

0.0, 

0.0, 

0.0 

3F+00, 

o.?e+oo. 

o 

. 

c 


O.O, 

O.O, 


0, 


0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 

O.O, 

0 . 0 , 

0.0, 

0.0, 

O.O, 


0, 

0.0, 

0.0, 

0.0, 

0.0, 

0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 

0 . 0 , 


0. lE^'Ol , 
0.6AE+00 
0.36F+00 
0. l6F*fr-00 
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QC = 4.62775674E+00 QCP = 4. 6067269?F^'00 OR = 0. QRP = 0. ORR = 9.89457607F-01 
HE == 3.72256809E+-0? HEP = 3. 71 712062E1-02 HWH == 3. 74851284E+02 HW ^ 3. 7341 8743F+02 X = 2 . 43 1 86 140E-02 
QAER0=-3,22535368E-02 QAER0P=-2. 1 1 5 13472E-02 TF =0. TB = 0. ORAT = 4. 68 18 18 1 8E-01 
PRAT = 1.60000000E-01 WTCI ^ 3 . 90000000F4- 00 ARFA = 2. OOOOOOOOF-01 WTBP == 4. 24776969E + 00 WTSliM= 4. 06104930E+00 
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TABLE I.- COMPARISON OF EXACT AND CALCULATED TEMPERATURES 
FOR CONSTANT PROPERTY SLABa 


Number 
of layers 

Error criterion 

(error < t' - t 

Maximum error 
pexact - Tcalc) x M 
V "^exact / 

1 

0.001 

0.5 

1 

.0001 

.005 

4 

.001 

.6 

4 

.0001 

.005 

b4 

.001 

.76 


.0001 

.64 


^The total thickness is the same for all cases. In the first two cases, 
41 stations were used. For the remaining cases, the first layer had 11 sta- 
tions and all other layers had 10 stations except that when the air-gap equa- 
tions were used, the third layer had only 1 station. 

t>The air-gap equations are used in the third layer. 
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Figure 2.- Simplified program flow chart. 
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Figure 3.- Computer program model of sample case heat shield system. 
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Insulation layer thickness, cm 

Figure 4.- Insulation layer thickness iterations for sample case. 





Percent vehicle area 


Figure 6.- Assumed heating rate distribution over vehicle. 
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Percent vehicle area 

Figure 7.- Cumulative average heat shield unit mass for 
assumed heating rate distribution. 
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